% this is  a general  rib waveguide 


clear all

import com.comsol.model.*
import com.comsol.model.util.*
global model;

model = ModelUtil.create(['Model' num2str(randi(200,1,1))]); 
model.name('rib_waveguide mode analysis');

geom1 = model.geom.create('geom1', 2);
model.geom('geom1').lengthUnit([native2unicode(hex2dec('00b5'), 'Cp1252') 'm']);  % set length unit in um

MyGeometry=geom1;DoSetMaterials=false;DoSetMesh=true;MyName00='A';
Geometry2D_Rib

MyNewMaterial=setMaterial(MyMaterial,'geom1',2 );
%%%%%%%


%  create mesh
global myMesh
myMesh=model.mesh.create('mesh1', 'geom1');
meshSize=myMesh.feature('size');
meshSize.set('hauto', '2');
%myMesh.feature.create('ftri1', 'FreeTri');
setMesh(MyMaterial,'geom1',2 );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.physics.create('emw', 'ElectromagneticWaves', 'geom1');

model.study.create('std1');
model.study('std1').feature.create('mode', 'ModeAnalysis');

model.physics('emw').feature('wee1').set('DisplacementFieldModel', 1, 'RefractiveIndex');
model.physics('emw').feature('wee1').set('n_mat', 1, 'from_mat');

model.study('std1').feature('mode').set('geomselection', 'geom1');
model.study('std1').feature('mode').set('physselection', 'emw');
model.study('std1').feature('mode').set('shift', 'nAprox');  %  aproximate mode refractive index
model.study('std1').feature('mode').set('modeFreq', ['c_const/' num2str(lambda) '[um]']);



model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature('st1').set('study', 'std1');
model.sol('sol1').feature('st1').set('studystep', 'mode');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature.create('e1', 'Eigenvalue');
model.sol('sol1').feature('e1').set('shift',  'nAprox');   %  aproximate mode refractive index
model.sol('sol1').feature('e1').set('neigs', numberModeSearch);                            %  number of mode to search
model.sol('sol1').feature('e1').set('transform', 'effective_mode_index');
model.sol('sol1').feature('e1').set('control', 'mode');
model.sol('sol1').feature('e1').feature('aDef').set('complexfun', true);
model.sol('sol1').attach('std1');

model.sol('sol1').runAll;
disp('calculations OK')

AreaAll=mphint(model,'1','Solnum','1');  % integration over all domains
AreaCore=mphint(model,'1','Solnum','1','selection',[MyNewMaterial(1).Selection; MyNewMaterial(2).Selection]);  %integrating over core and clad
% find which solutions are realistic
n_eff_all = mphglobal(model,'emw.neff','Complexout','on');
jMode=0;
for j=1:length(n_eff_all)
    %if real(n_eff_all(j))>nAlGaAs(xClad,lambda)      
        jMode=jMode+1;
        mode(jMode).neff=real(n_eff_all(j));
        mode(jMode).loss=kdb(imag(n_eff_all(j)),lambda)/10;  %loss in dB/mm
       
        mode(jMode).n=n_eff_all(j);
        mode(jMode).solutionN=j;
   % end
end
disp([' mode total number:  ' num2str(jMode)])



%%%%%% find whether the mode TE or TM
jTM=0;jTE=0;jTEM=0;

for j=1:jMode
    
mode(j).TEintens=mphint(model,'abs(emw.Ex)^2','Solnum',num2str(mode(j).solutionN));  % integration over all domains
mode(j).TMintens=mphint(model,'abs(emw.Ey)^2','Solnum',num2str(mode(j).solutionN));
if mode(j).TEintens>3*mode(j).TMintens 
    mode(j).Polar='TE';
    jTE=jTE+1;TEmod(jTE,1)=mode(j).solutionN;TEmod(jTE,2)=mode(j).neff;
    
elseif mode(j).TMintens>3*mode(j).TEintens
    mode(j).Polar='TM';
    jTM=jTM+1;TMmod(jTM,1)=mode(j).solutionN;TMmod(jTM,2)=mode(j).neff;
else
    mode(j).Polar='TEM';
    jTEM=jTEM+1;TEMmod(jTEM,1)=mode(j).solutionN;TEMmod(jTEM,2)=mode(j).neff;
end

end

% put each mode in order, mode with highest refractive index first
mes='';
if jTE>0 TEmod=sortrows(TEmod,-2); mes=['number of TE modes ' num2str(size(TEmod,1))];end
if jTE>0 TMmod=sortrows(TMmod,-2);mes=[mes ' of TM modes ' num2str(size(TMmod,1))]; end
 if jTEM>0 TEMmod=sortrows(TEMmod,-2);mes=[mes ' of TEM modes ' num2str(size(TEMmod,1))]; end

disp(mes )

%%%%%%%%%%%% plot result
jplot=3;
  
if jTE>0
for j=1:size(TEmod,1)
    jplot=jplot+1;
dataName=['pg' num2str(jplot)];
model.result.create(dataName, 'PlotGroup2D');
model.result(dataName).feature.create('surf1', 'Surface');
model.result(dataName).feature('surf1').set('descr', ['Electric field TE' num2str(j-1) ' mode, x component']);
model.result(dataName).feature('surf1').set('expr', 'abs(emw.Ex)');
model.result(dataName).name(['TE' num2str(j-1)]);
model.result(dataName).set('solnum', num2str(TEmod(j,1)));
%model.result('dataName').set('data', 'dset1');
model.result(dataName).run;
end 
end

if jTM>0
for j=1:size(TMmod,1)
    jplot=jplot+1;
dataName=['pg' num2str(jplot)];
model.result.create(dataName, 'PlotGroup2D');
model.result(dataName).feature.create('surf1', 'Surface');
model.result(dataName).feature('surf1').set('descr', ['Electric field TE' num2str(j-1) ' mode, x component']);
model.result(dataName).feature('surf1').set('expr', 'abs(emw.Ey)');
model.result(dataName).name(['TM' num2str(j-1)]);
model.result(dataName).set('solnum', num2str(TMmod(j,1)));
%model.result(dataName).set('data', 'dset1');
model.result(dataName).run;
end 
end



disp('all OK')